! $Id: test_any_kernel.F90,v 1.7 2008/11/07 20:50:34 nnz Exp nnz $
!
!This is a test driver to validate the kernels
!In order to compile, first make the kernel libraries upstairs (make -C ../ ), then do
!
! ifort -g -r8 -i4 test_any_kernel.F90 ../libHOME_kernels.a -I../ -o test_any_kernel
!
! gfortran test_any_kernel.F90 ../libHOME_kernels.a -I../ -o test_any_kernel 
!
!To run: ./test_any_kernel
!To view the advected profiles:  xmgrace xTracerFlx_mdflSsupB_cne.dat
!
!The Test Box Boundaries extends from -NCP to +NCP (so 2*NCP+1 points) with NHALO halos at each side

#define NCP   25
#define NHALO 4
!Number of time intervals
#define TN 410 !410    

!----------------------------
!Data and Compute indices
!
#define IND  2*NCP+1+2*NHALO   
#define ISD -NCP-NHALO     
#define IED  NCP+NHALO   
#define ISC -NCP     
#define IEC  NCP   

#define JND  2*NCP+1+2*NHALO   
#define JSD -NCP-NHALO   
#define JED  NCP+NHALO   
#define JSC -NCP     
#define JEC  NCP   

!shape_vector dimmensions
#define SV1 2
#define SV2 8

program test_kernels

  use HOME_kernels_cne
  use HOME_kernels_csw


  implicit none

  integer :: i,j,k,nt
  real    :: dt, error, maxError,xt
  real    :: period,rms
  real,parameter    :: tol=1.0e-11
  integer :: shape_vector(SV1,SV2) 
  real   , dimension(ISD:IED,JSD:JED) :: u,v,uhrho,vhrho,dx,dy,tmask, tracer_flux
  real   , dimension(ISD:IED,JSD:JED,0:TN,-5:5) :: tracer
  integer :: nprofiles,p,p0,ker, error_flag
  character(len=80),dimension(62) :: kernel_name = (/&
       trim('xTracerFlx_pqm0_cne'),&
       trim('xTracerFlx_ppm0_cne'),&
       trim('yTracerFlx_ppm0_cne'),&
       trim('xTracerFlx_ppm0_csw'),&
       trim('yTracerFlx_ppm0_csw'),&
       trim('yTracerFlx_pqm0_cne'),&
       trim('xTracerFlx_pqm0_csw'),&
       trim('yTracerFlx_pqm0_csw'),&
       trim('xTracerFlx_pqm1_cne'),&
       trim('yTracerFlx_pqm1_cne'),&
       trim('xTracerFlx_pqm1_csw'),&
       trim('yTracerFlx_pqm1_csw'),&
       trim('xTracerFlx_upwind_cne'),& 
       trim('yTracerFlx_upwind_cne'),& 
       trim('xTracerFlx_upwind_csw'),& 
       trim('yTracerFlx_upwind_csw'),& 
       trim('xTracerFlx_mdflSupB_cne'),&
       trim('yTracerFlx_mdflSupB_cne'),&
       trim('xTracerFlx_mdflSupB_csw'),&
       trim('yTracerFlx_mdflSupB_csw'),&
       trim('xTracerFlx_plm_cne'),&
       trim('yTracerFlx_plm_cne'),&
       trim('xTracerFlx_plm_csw'),&
       trim('yTracerFlx_plm_csw'),&
       trim('xTracerFlx_ppm1_cne'),&
       trim('yTracerFlx_ppm1_cne'),&
       trim('xTracerFlx_ppm1_csw'),&
       trim('yTracerFlx_ppm1_csw'),&
       trim('xTracerFlx_ppm2_cne'),&
       trim('yTracerFlx_ppm2_cne'),&
       trim('xTracerFlx_ppm2_csw'),&
       trim('yTracerFlx_ppm2_csw'),&
       trim('xTracerFlx_ppm3_cne'),&
       trim('yTracerFlx_ppm3_cne'),&
       trim('xTracerFlx_ppm3_csw'),&
       trim('yTracerFlx_ppm3_csw'),&
       trim('xTracerFlx_DST3FL_cne'), &
       trim('xTracerFlx_DST3FL_csw'), &
       trim('yTracerFlx_DST3FL_cne'), &
       trim('yTracerFlx_DST3FL_csw'), &
       trim('xTracerFlx_OS7MP_cne'), &
       trim('xTracerFlx_OS7MP_csw'), &
       trim('yTracerFlx_OS7MP_cne'), &
       trim('yTracerFlx_OS7MP_csw'), &
       trim('xTracerFlx_plmGOLD_cne'),&
       trim('yTracerFlx_plmGOLD_cne'),&
       trim('xTracerFlx_plmGOLD_csw'),&
       trim('yTracerFlx_plmGOLD_csw'),&
       trim('xTracerFlx_laxwen_cne'), &
       trim('xTracerFlx_laxwen_csw'), &
       trim('yTracerFlx_laxwen_cne'), &
       trim('yTracerFlx_laxwen_csw'), &
       trim('xTracerFlx_plmLaurent_cne'), &
       trim('xTracerFlx_plmLaurent_csw'), &
       trim('xTracerFlx_2ndOrderCentered_cne'),& 
       trim('yTracerFlx_2ndOrderCentered_cne'),& 
       trim('xTracerFlx_2ndOrderCentered_csw'),& !has x--> -x problem
       trim('yTracerFlx_2ndOrderCentered_csw'),& !has y--> -y problem
       trim('xTracerFlx_4thOrderCentered_cne'),&
       trim('yTracerFlx_4thOrderCentered_cne'),&
       trim('xTracerFlx_4thOrderCentered_csw'),&
       trim('yTracerFlx_4thOrderCentered_csw')&
       /)

  !Error file
  open(16,file='error_out',status='unknown')
  open(41,file='report',status='unknown')

  shape_vector(1,1) = ISD
  shape_vector(1,2) = IED
  shape_vector(1,3) = ISD
  shape_vector(1,4) = IED
  shape_vector(1,5) = ISD
  shape_vector(1,6) = IED

  shape_vector(2,1) = JSD
  shape_vector(2,2) = JED
  shape_vector(2,3) = JSD
  shape_vector(2,4) = JED
  shape_vector(2,5) = JSD
  shape_vector(2,6) = JED


  !Velocities and grid sizes
  do i=ISD,IED ; do j=JSD,JED
     u(i,j) = 1
     v(i,j) = 1
     uhrho(i,j) = 1
     vhrho(i,j) = 1
     dx(i,j) = 1.
     dy(i,j) = 1.
     tmask(i,j) = 1.0
  enddo; enddo

  !Island

  !Some of the kernels require properly haloed tmask
  call halosync(tmask)

  !Notes: 
  !Forward time and one-sided space scheme is used to obtain the difference equation here.
  !For u=1 the numerical solution amplitude is constant and produces the true solution.
  !For u<1 the numerical solution amplitude is damping.
  !For u>1 the numerical solution amplitude is growing.

  dt=0.25 ! ==> with the above choices for u we have cfl=dt

  period=(2*NCP+1)/dt  ! Why +1? It should be period= L/cfl = 2*NCP/dt 


  do ker=1,2 ! 62

     error_flag=0

     write(6,*)  'Testing advection kernel ', trim(kernel_name(ker))
     write(16,*)  'Testing advection kernel ', trim(kernel_name(ker))
     write(16,*) '         ker           p           nt  maxError'

     !
     !Initial tracer concentration profile
     !
     call create_initial_profiles(tracer,nprofiles)

     if(INDEX(kernel_name(ker), '_cne').ne.0) then
        shape_vector(1,7) = ISC-1 ! -1 needed for cne
        shape_vector(1,8) = IEC   ! +1 needed for csw
        shape_vector(2,7) = JSC-1 ! -1 needed for cne
        shape_vector(2,8) = JEC   ! +1 needed for csw
     elseif(INDEX(kernel_name(ker), '_csw').ne.0) then
        shape_vector(1,7) = ISC   ! -1 needed for cne
        shape_vector(1,8) = IEC+1 ! +1 needed for csw
        shape_vector(2,7) = JSC   ! -1 needed for cne
        shape_vector(2,8) = JEC+1 ! +1 needed for csw
     endif
     !
     !Advect each profile independently
     !
     do p=-nprofiles,nprofiles

        do nt=0,TN-1
           !zero out the tracer_flux
           tracer_flux = 0.0       
           !call a kernel to calculate tracer_flux
           !For later comparison advect positive index profiles in + direction and negative index profiles in - direction
           if(p.ge.0) then
              call some_kernel(kernel_name(ker),shape_vector,tmask,dt,dx,dy, u, v, uhrho, vhrho,tracer(:,:,nt,p),&
                   tracer_flux)
           else
              call some_kernel(kernel_name(ker),shape_vector,tmask,dt,dx,dy,-u,-v,-uhrho,-vhrho,tracer(:,:,nt,p),&
                   tracer_flux)
           endif

           !simple forward time step each of the profiles
           if(INDEX(kernel_name(ker), 'x').eq.1) then
              if(INDEX(kernel_name(ker), '_cne').ne.0) then
                 do j=JSC,JEC ; do i=ISC,IEC  
                    tracer(i,j,nt+1,p) = tracer(i,j,nt,p) - (tracer_flux(i,j) - tracer_flux(i-1,j)) * dt &
                         / dx(i,j) / dy(i,j) 
                 enddo; enddo
              elseif(INDEX(kernel_name(ker), '_csw').ne.0) then
                 do j=JSC,JEC ; do i=ISC,IEC  
                    tracer(i,j,nt+1,p) = tracer(i,j,nt,p) - (tracer_flux(i+1,j) - tracer_flux(i,j)) * dt &
                         / dx(i,j) / dy(i,j) 
                 enddo; enddo
              else
                 write(6,*) 'Error: Unknown kernel class ', kernel_name(ker)
              endif
           elseif(INDEX(kernel_name(ker), 'y').eq.1) then
              if(INDEX(kernel_name(ker), '_cne').ne.0) then
                 do j=JSC,JEC ; do i=ISC,IEC  
                    tracer(i,j,nt+1,p) = tracer(i,j,nt,p) - (tracer_flux(i,j) - tracer_flux(i,j-1)) * dt &
                         / dx(i,j) / dy(i,j) 
                 enddo; enddo
              elseif(INDEX(kernel_name(ker), '_csw').ne.0) then
                 do j=JSC,JEC ; do i=ISC,IEC  
                    tracer(i,j,nt+1,p) = tracer(i,j,nt,p) - (tracer_flux(i,j+1) - tracer_flux(i,j)) * dt &
                         / dx(i,j) / dy(i,j) 
                 enddo; enddo
              else
                 write(6,*) 'Error: Unknown kernel class ', kernel_name(ker)
              endif
           else
              write(6,*) 'Error: Unknown kernel direction ', kernel_name(ker)
           endif


        enddo
     enddo

     !
     ! Tests
     !
     ! If u ---> -u .AND. x --> A-x advection equation remains the same
     ! Hence if we prepare two initial profiles t_p and t_m such that t_m(x) = t_p(A-x) 
     ! and let them propagate in +u and -u directions respectively then we should get
     ! t_m(i,j;t) = t_p(A-i,j;t) for each t 
     ! We choose A=0 for simplicity and compare the positive index profiles with negative index ones
     ! Dito x--> y and u--> v
     !
     do p=1,nprofiles
        do nt=1,TN-1
           maxError = 0
           if(INDEX(kernel_name(ker), 'x').eq.1) then
              j=0
              do i=ISC+1,IEC-1     
                 error=abs(tracer(i,j,nt,p)-tracer(-i,j,nt,-p))
                 if(error.gt.maxError) maxError = error
              enddo
           elseif(INDEX(kernel_name(ker), 'y').eq.1) then
              i=0
              do j=JSC+1,JEC-1     
                 error=abs(tracer(i,j,nt,p)-tracer(i,-j,nt,-p))
                 if(error.gt.maxError) maxError = error
              enddo
           endif

           write(16,*) ker,p,nt, maxError
           if(maxError.gt.tol) error_flag=1
        enddo
     enddo

     !
     !Plots
     !
     open(40,file=trim(kernel_name(ker))//'.dat',status='unknown')
     do p=0,nprofiles
        do nt=0,TN,period
           if(INDEX(kernel_name(ker), 'x').eq.1) then
              j=0
              do i=ISC,IEC  
                 write(40,*) i+60*p,tracer(i,j,nt,p)
              enddo
              write(40,*) '&'
           elseif(INDEX(kernel_name(ker), 'y').eq.1) then
              i=0
              do j=JSC,JEC  
                 write(40,*) j,tracer(i,j,nt,p)
              enddo
              write(40,*) '&'
           endif
        enddo
     enddo
     close(40)

     !Calculate RMS at each period
     write(41,*) '    profile      period  rms'
     write(41,*) trim(kernel_name(ker))
     do p=0,nprofiles
        do nt=period,TN,period
           rms=0
           do i=ISC,IEC  
              rms=rms+(tracer(i,0,nt,p) - tracer(i,0,0,p))**2
           enddo
           rms = sqrt(rms)
           write(41,*) p,int(nt/period),rms
        enddo
     enddo

     if(error_flag.ne.0) then
        write(6,*) trim(kernel_name(ker)), ' FAILED with error flag ', error_flag
        write(41,*) trim(kernel_name(ker)), ' FAILED with error flag ', error_flag
     else
        write(6,*) trim(kernel_name(ker)), ' Passed'
        write(41,*) trim(kernel_name(ker)), ' Passed'
     endif
     write(6,*) ''
     write(41,*) ''

     !Exact Advection ?

     !  j=0
     !  do p=2,2 !-nprofiles,nprofiles
     !     do nt=0,TN-1,TN/4
     !        xt=nt*dt
     !        do i=ISC,IEC  
     !           write(51,*) i-xt,tracer(i,j,nt,p)
     !        enddo        
     !        write(51,*) '&'
     !        do i=ISC,IEC  
     !           write(52,*) i-xt,  exp(-(i-xt+j)**2 / 8.)
     !        enddo        
     !        write(52,*) '&'
     !     enddo
     !  enddo

     !Compare advected profiles
     if(trim(kernel_name(ker)) .eq. 'xTracerFlx_pqm1_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_ppm3_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_ppm2_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_ppm1_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_ppm0_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_pqm0_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_ppmNiki_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_laxwen_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_plmGOLD_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_plmLaurent_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_OS7MP_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_DST3FL_cne' .or. &
          trim(kernel_name(ker)) .eq. 'xTracerFlx_plm_cne') then  

        j=0
        nt=period
        do p=0,nprofiles           
           write(50+p,*) '#',trim(kernel_name(ker))
           do i=ISC,IEC  
              write(50+p,*) i,tracer(i,j,nt,p)
           enddo
           write(50+p,*) '&'        
        enddo
     endif

  enddo !do ker=1,16


  do p=0,nprofiles           
     write(50+p,*) '#Exact Profile'
     do i=ISC,IEC  
        write(50+p,*) i,tracer(i,0,0,p)
     enddo
  enddo

contains

  subroutine some_kernel(kernel_name, shape_vector,tmask,dt,dx,dy,u,v,uhrho,vhrho,tracer_field,tracer_flux)
    character(len=*), intent(in) :: kernel_name
    integer,          intent(in) :: shape_vector(SV1,SV2)       ! Shape vector
    real,             intent(in) :: tmask(ISD:IED,JSD:JED)    	! tracer mask
    real,             intent(in) :: dt 
    real,             intent(in) :: dx(ISD:IED,JSD:JED)    	! x width of face at u-point (m)
    real,             intent(in) :: dy(ISD:IED,JSD:JED)    	! y width of face at u-point (m)
    real,             intent(in) :: u(ISD:IED,JSD:JED)  	! x Component of velocity normal to face (Kg/m/s)
    real,             intent(in) :: v(ISD:IED,JSD:JED)  	! y Component of velocity normal to face (Kg/m/s)
    real,             intent(in) :: uhrho(ISD:IED,JSD:JED)  	! x Component of mass flow normal to face (Kg/m/s)
    real,             intent(in) :: vhrho(ISD:IED,JSD:JED)  	! y Component of mass flow normal to face (Kg/m/s)
    real,             intent(inout) :: tracer_field(ISD:IED,JSD:JED)	! Tracer field at c-points ([tracer_concentration]) 
    real,             intent(out):: tracer_flux(ISD:IED,JSD:JED)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)


    !Some of the kernels require properly haloed tracer
    call halosync(tracer_field)


    select case(kernel_name)
    case ('xTracerFlx_upwind_cne') 
       call xTracerFlx_upwind_cne(shape_vector, dy, uhrho,tracer_field,tracer_flux)     
    case ('yTracerFlx_upwind_cne') 
       call yTracerFlx_upwind_cne(shape_vector, dx, vhrho,tracer_field,tracer_flux)     
    case ('xTracerFlx_upwind_csw') 
       call xTracerFlx_upwind_csw(shape_vector, dy, uhrho,tracer_field,tracer_flux)     
    case ('yTracerFlx_upwind_csw') 
       call yTracerFlx_upwind_csw(shape_vector, dx, vhrho,tracer_field,tracer_flux)     
    case ('xTracerFlx_2ndOrderCentered_cne') 
       call xTracerFlx_2ndOrderCentered_cne(shape_vector,dy,uhrho,tracer_field,tracer_flux)
    case ('yTracerFlx_2ndOrderCentered_cne') 
       call yTracerFlx_2ndOrderCentered_cne(shape_vector,dx,vhrho,tracer_field,tracer_flux)
    case ('xTracerFlx_2ndOrderCentered_csw') 
       call xTracerFlx_2ndOrderCentered_csw(shape_vector,dy,uhrho,tracer_field,tracer_flux)
    case ('yTracerFlx_2ndOrderCentered_csw') 
       call yTracerFlx_2ndOrderCentered_csw(shape_vector,dx,vhrho,tracer_field,tracer_flux)
    case ('xTracerFlx_4thOrderCentered_cne') 
       call xTracerFlx_4thOrderCentered_cne(shape_vector,tmask,dy,uhrho,tracer_field,tracer_flux)
    case ('yTracerFlx_4thOrderCentered_cne') 
       call yTracerFlx_4thOrderCentered_cne(shape_vector,tmask,dx,vhrho,tracer_field,tracer_flux)
    case ('xTracerFlx_4thOrderCentered_csw') 
       call xTracerFlx_4thOrderCentered_csw(shape_vector,tmask,dy,uhrho,tracer_field,tracer_flux)
    case ('yTracerFlx_4thOrderCentered_csw') 
       call yTracerFlx_4thOrderCentered_csw(shape_vector,tmask,dx,vhrho,tracer_field,tracer_flux)
    case ('xTracerFlx_mdflSupB_cne') 
       call xTracerFlx_mdflSupB_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_mdflSupB_cne') 
       call yTracerFlx_mdflSupB_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_mdflSupB_csw') 
       call xTracerFlx_mdflSupB_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_mdflSupB_csw') 
       call yTracerFlx_mdflSupB_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plm_cne') 
       call xTracerFlx_plm_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_plm_cne') 
       call yTracerFlx_plm_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plm_csw') 
       call xTracerFlx_plm_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_plm_csw') 
       call yTracerFlx_plm_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm0_cne') 
       call xTracerFlx_ppm_cne(shape_vector,0,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm0_cne') 
       call yTracerFlx_ppm_cne(shape_vector,0,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm0_csw') 
       call xTracerFlx_ppm_csw(shape_vector,0,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm0_csw') 
       call yTracerFlx_ppm_csw(shape_vector,0,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm1_cne') 
       call xTracerFlx_ppm_cne(shape_vector,1,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm1_cne') 
       call yTracerFlx_ppm_cne(shape_vector,1,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm1_csw') 
       call xTracerFlx_ppm_csw(shape_vector,1,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm1_csw') 
       call yTracerFlx_ppm_csw(shape_vector,1,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm2_cne') 
       call xTracerFlx_ppm_cne(shape_vector,2,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm2_cne') 
       call yTracerFlx_ppm_cne(shape_vector,2,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm2_csw') 
       call xTracerFlx_ppm_csw(shape_vector,2,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm2_csw') 
       call yTracerFlx_ppm_csw(shape_vector,2,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm3_cne') 
       call xTracerFlx_ppm_cne(shape_vector,3,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm3_cne') 
       call yTracerFlx_ppm_cne(shape_vector,3,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppm3_csw') 
       call xTracerFlx_ppm_csw(shape_vector,3,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_ppm3_csw') 
       call yTracerFlx_ppm_csw(shape_vector,3,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plmLaurent_cne') 
       call xTracerFlx_plmLaurent_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plmLaurent_csw') 
       call xTracerFlx_plmLaurent_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_OS7MP_cne') 
       call xTracerFlx_OS7MP_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_OS7MP_cne') 
       call yTracerFlx_OS7MP_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_OS7MP_csw') 
       call xTracerFlx_OS7MP_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_OS7MP_csw') 
       call yTracerFlx_OS7MP_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_DST3FL_cne') 
       call xTracerFlx_DST3FL_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_DST3FL_cne') 
       call yTracerFlx_DST3FL_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_DST3FL_csw') 
       call xTracerFlx_DST3FL_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_DST3FL_csw') 
       call yTracerFlx_DST3FL_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plmGOLD_cne') 
       call xTracerFlx_plmGOLD_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_plmGOLD_cne') 
       call yTracerFlx_plmGOLD_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_plmGOLD_csw') 
       call xTracerFlx_plmGOLD_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_plmGOLD_csw') 
       call yTracerFlx_plmGOLD_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_ppmNiki_cne') 
       call xTracerFlx_ppmNiki_cne(shape_vector,0,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_pqm0_cne') 
       call xTracerFlx_pqm_cne(shape_vector,0,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_pqm1_cne') 
       call xTracerFlx_pqm_cne(shape_vector,1,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_pqm0_cne') 
       call yTracerFlx_pqm_cne(shape_vector,0,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_pqm1_cne') 
       call yTracerFlx_pqm_cne(shape_vector,1,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_pqm0_csw') 
       call xTracerFlx_pqm_csw(shape_vector,0,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_pqm1_csw') 
       call xTracerFlx_pqm_csw(shape_vector,1,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_pqm0_csw') 
       call yTracerFlx_pqm_csw(shape_vector,0,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_pqm1_csw') 
       call yTracerFlx_pqm_csw(shape_vector,1,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_laxwen_cne') 
       call xTracerFlx_laxwen_cne(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_laxwen_cne') 
       call yTracerFlx_laxwen_cne(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case ('xTracerFlx_laxwen_csw') 
       call xTracerFlx_laxwen_csw(shape_vector,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)
    case ('yTracerFlx_laxwen_csw') 
       call yTracerFlx_laxwen_csw(shape_vector,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)
    case default
       write(6,*) 'Error: The specified kernel ',trim(kernel_name),' is not implemented!'
       stop
    end select

  end subroutine some_kernel

  subroutine halosync0(q)
    ! Arguments
    real, intent(inout) :: q(ISD:IED,JSD:JED)
    ! Local variables
    integer :: h

    do h=1,NHALO
       q(ISC-h,:)=q(IEC-h,:)
       q(IEC+h,:)=q(ISC+h,:)
    enddo

    do h=1,NHALO
       q(:,JSC-h)=q(:,JEC-h+1)
       q(:,JEC+h)=q(:,JSC+h-1)
    enddo

  end subroutine halosync0

  subroutine halosync(q)
    ! Arguments
    real, intent(inout) :: q(ISD:IED,JSD:JED)
    ! Local variables
    integer :: h

    do h=1,NHALO
       q(ISC-h,:)=q(IEC-h+1,:)
       q(IEC+h,:)=q(ISC+h-1,:)
    enddo

    do h=1,NHALO
       q(:,JSC-h)=q(:,JEC-h+1)
       q(:,JEC+h)=q(:,JSC+h-1)
    enddo

  end subroutine halosync

  subroutine create_initial_profiles(tracer,nprofiles)
    real   , dimension(ISD:IED,JSD:JED,0:TN,-5:5), intent(inout) :: tracer
    integer, intent(out) :: nprofiles
    integer :: i,j,p

    tracer = 0.0

    !double peaks on p=0
    p=0

    do j=-10,10; do i=-10,10
       tracer(i,j,0,p) = exp(-(i+j)**2 / 8.)  
    enddo; enddo

    do j=JSD,JED ; do i=15,20
       tracer(i,j,0,p) = 1.0/2 
    enddo; enddo
    do j=15,20 ; do i=ISD,IED
       tracer(i,j,0,p) = 1.0/2 
    enddo; enddo

    !step for p=1

    do i=-10,10; do j=-10,10  
       tracer(i,j,0,1) = 1.0 
    enddo; enddo

    !Smooth peak on p=2,-2
    do i=-20,20; do j=-20,20  
       tracer(i,j,0,2) = exp(-(i+j)**2 / 15.) 
    enddo; enddo

    !Smooth peak on p=3,-3
    ! do i=-9,9; do j=-9,9  
    !    tracer(i,j,0,3) = 1-((i+j)/10.)**4 
    ! enddo; enddo

    do j=JSD,JED ; do i=ISD,IED  
       tracer(i,j,0,3) = (((i+j)/10.)**2 - 1)*(((i+j)/20.)**2 - 1) 
    enddo; enddo



    nprofiles = 3
    do p=1,nprofiles
       tracer(:,:,0,-p) = tracer(:,:,0,p) 
    enddo

    !nprofiles = 0
  end subroutine create_initial_profiles

end program test_kernels
